function cal (L)

   ## usage: cal (L)
   ## 
   ## calibrate model

   global PAR ADIR PRJ clim eof

   pkg load general io nan

   ldata = fullfile("data", PRJ, PAR.pdd) ;

   ofile = fullfile(ldata, "pdd.ob") ;
   vfile = fullfile(ADIR, [strrep(PAR.ana, "/", "_") ".pc.ob"]) ;
   if (isnewer(PAR.mfile, ofile, vfile) && !L), return ; endif

   f = dbstack ; printf("\n\t###   %s   ###\n\n", f(1).name) ;

   ws = sprintf("    Calibration of XDS using %s    ", PAR.sim{1}) ; nws = length(ws) ;
   wt1 = "#"(ones(1,7)) ; wt2 = "#"(ones(1,nws)) ;
   printf(["\n\t" wt1 wt2 wt1 "\n"]) ;
   printf(["\t" wt1 ws wt1 "\n"]) ;
   printf(["\t" wt1 wt2 wt1 "\n\n"]) ;

   ## get clim and pdd
   load(cfile = fullfile(ldata, "clim.ob")) ;
   printf("cal: <-- %s\n", cfile) ;
   load(ofile) ;
   printf("cal: <-- %s\n", ofile) ;
   pdd = do_cov(pdd) ;

   ## get pc of ptr, define X and Y
   load(vfile) ;
   printf("cal: <-- %s\n", vfile) ;

   X = selper(PC, sdate(PC.id, PAR.cper(1,:), PAR.cper(2,:))) ;
   Y = selper(pdd, sdate(pdd.id, PAR.cper(1,:), PAR.cper(2,:))) ;
   [X.nr X.nc] = size(X.x) ;
   X.CI = sdate(X.id, PAR.xdscal(1,:), PAR.xdscal(2,:)) ;
   if isfield(PAR, "xdsval")
      X.VI = sdate(X.id, PAR.xdsval(1,:), PAR.xdsval(2,:)) ;
   else
      X.VI = !X.CI ;
   endif
   if strcmp(PAR.optptr, "manual") && (!any(X.CI) || !any(X.VI))
      error("cal: too few values for cal/val") ;
   endif 

   if isfield(PAR, "ptridx") && !isempty(PAR.ptridx)
      Jptr = false ;
      for j=1:length(PAR.ptridx)
	 Jptr = Jptr | X.J == PAR.ptridx(j) ;
      endfor 
   else
      Jptr = true(1, X.nc) ;
   endif
   Xj = Xjfun(X.J) ;

   if PAR.pcr
      X.x(:,!Jptr) = nan ;
      ## [eof.E eof.PC eof.ev] = pca(X.x, PAR.trunc) ;
      [eof.E eof.PC eof.ev] = pca(X.x) ;
      eof.E = nan2lim(eof.E, 0) ;
      X.x = X.x * eof.E ;
      [X.nr X.nc] = size(X.x) ; Jptr = true(1, X.nc) ;
      printf("%s: --> %s\n", f(1).name, efile = strrep(PAR.mfile, "mod", "pcr")) ; 
      save(efile, "eof") ;
   endif

   ## select predictors
   printf("%s: Using %d potential predictors (PCs).\n", f(1).name, sum(Jptr)) ;
   Jeof = ptrsel(X, Y, Jptr) ;
   fid = fopen(PAR.plist, "w") ; J = find(Jeof) ;
   if PAR.pcr
      for j=1:sum(Jeof)
	 fprintf(fid, "%d: EOF %d\n", j, J(j)) ;
      endfor
   else
      for j=1:sum(Jeof)
	 fprintf(fid, "%d: %s %d\n", j, X.ptr(X.J){J(j)}, Xj(J(j))) ;
      endfor
   endif
   fclose(fid) ;

   skl(Jeof, Y, X, "disp") ;

   ## best model \\
   X = selper(PC, sdate(PC.id, PAR.cper(1,:), PAR.cper(2,:))) ;
   if PAR.pcr
      X.x = X.x * eof.E ;
      [X.nr X.nc] = size(X.x) ;   
   endif
   [E R] = do_cal(Y, X, Jeof, PAR.pcr) ;
   printf("cal: --> %s\n", PAR.mfile) ; 
   save(PAR.mfile, "E", "R", "Jeof") 
   ## best model //

endfunction


function res = skl (Jeof, Y, X, mode)

   ## usage:  res = skl (Jeof, Y, X, mode)
   ##
   ## validate vars rel. to number of EOFs

   global PAR clim

   wX = X ; wY = Y ; wX.x(!X.CI,:) = nan ; wY.z(!X.CI,:) = nan ;
   [E R] = do_cal(wY, wX, Jeof) ;

   ## validation
   y = nan(rows(X.x), size(E, 2)) ;
   for is = 1:length(PAR.ssn)
      I = selssn(PAR.ssn{is}, X.id) ;
      y(I,:) = nanmult(X.x(I,:), E(:,:,is)) ;
   endfor

   if !isempty(PAR.pddadj) && strcmp(PAR.mod, "xds")
      ncpu = PAR.ncpu ; # ugly
      PAR.ncpu = 1 ;
      PAR.dbg = false ;
      pp = feval(PAR.adjfun, y) ;
      y = feval(PAR.adjfun, y, pp) ;
      PAR.ncpu = ncpu ;
   endif

   VAR = unique(Y.vars) ; cres = vres = [] ; dVAR = {} ;
   J0 = ind2log(PAR.validx, Y.nc) ;

   for v = VAR
      J = strcmp(Y.vars, v{:}) & J0 ;
      if !any(J), continue ; endif
      dVAR = union(dVAR, v) ;
      [w(1) w(2) w(3) w(4)] = re(mean(y(X.CI,J), 2), mean(Y.z(X.CI,J), 2)) ; cres = [cres w(PAR.valscore)] ;
      [w(1) w(2) w(3) w(4)] = re(mean(y(X.VI,J), 2), mean(Y.z(X.VI,J), 2)) ; vres = [vres w(PAR.valscore)] ;
   endfor
   res = [cres ; vres] ;
   if nargin > 3
      f = dbstack ;
      score = {"EV" "correlation" "mean bias" "std bias"} ;
      printf("%s: calibration/validation %s skill ([%%]):\n", f(1).name, score{PAR.valscore}) ;
      displaytable(res', {"cal" "val"}, 5, {".0f"}, dVAR, 1, " ") ;
   endif

   if !PAR.kstest, return ; endif 

   wrn = warning ;
   warning("off") ;

   cres = vres = [] ;
   for v = VAR
      J = strcmp(Y.vars, v{:}) & J0 ;
      if !any(J), continue ; endif
      wy = zscore(y(X.CI,J)) ;
      for j = 1:columns(wy)
	 ks(j) = kolmogorov_smirnov_test(wy(:,j), "norm", 0, 1) ;
      endfor 
      cres = [cres mean(ks)] ;
      wy = zscore(y(X.VI,J)) ;
      for j = 1:columns(wy)
	 ks(j) = kolmogorov_smirnov_test(wy(:,j), "norm", 0, 1) ;
      endfor 
      vres = [vres mean(ks)] ;
   endfor
   res = [res ; [cres ; vres]] ;

   warning(wrn) ;

endfunction


function neof = neof_fig (NEOF, res, VAR)

   ## usage:  neof = neof_fig (NEOF, res, VAR)
   ##
   ## estimate neof from figures

   global PAR

   if (PAR.valscore == 1)
      ylbl = "EV [%]" ;
   else
      ylbl = "rho" ;
   endif

   nc = 2 + 2*PAR.kstest ; nV = length(VAR) ;
   set(0, "defaulttextfontsize", 8, "defaultaxesfontsize", 8) ;

   for j = 1:nV
      v = VAR{j} ;
      subplot(nc,nV,j) ;
      plot(NEOF', squeeze(res(1,j,:))) ; set(gca, "title", ["cal(" v ")"], "xlabel", "neof", "ylabel", ylbl) ;
      subplot(nc,nV,nV+j) ;
      plot(NEOF', squeeze(res(2,j,:))) ; set(gca, "title", ["val(" v ")"], "xlabel", "neof", "ylabel", ylbl) ;
      if PAR.kstest
	 subplot(nc,nV,2*nV+j) ;
	 plot(NEOF', squeeze(res(3,j,:))) ; title(["KS(cal " v ")"]) ;
	 subplot(nc,nV,3*nV+j) ;
	 plot(NEOF', squeeze(res(4,j,:))) ; title(["KS(val " v ")"]) ;
      endif
   endfor
   if PAR.dbg, print([fileparts(PAR.mfile) "/cal.png"]) ; endif

   if !strcmp(getenv("TERM"), "dumb")
      printf("cal: How many predictors (shown on x-axis)?\n") ;
      printf("cal: Suggested: 3-10 x number of local variables,\n") ;
      printf("cal: use neof with %d <= NEOF <= %d\n", NEOF([1 end])) ;
      printf("cal: use \"<var>\" (with quotes!) for optimum var\n") ;

      in = input("cal: NEOF = ") ;
   else
      in = NEOF(end) ;
      sleep(3) ;
   endif
   close ;

   if isnumeric(in) && in >= NEOF(1) && in <= NEOF(end)
      neof = in ;
   elseif ischar(in)
      jv = find(strcmp(VAR, in)) ;      
      [w j] = max(res(2,jv,:)) ;
      neof = NEOF(j) ;
   else
      error("cal: NEOF invalid, exiting.\n")
   endif 

endfunction


function Jeof = ptrsel (X, Y, Jptr)

   ## usage: Jeof = ptrsel (X, Y, Jptr)
   ## 
   ## predictor selection

   global PAR
   f = dbstack ;

   VAR = unique(Y.vars) ;

   switch PAR.optptr

      case "auto"
	 I = !any(isnan([X.x Y.z]), 2) ; Jeof = false(1,X.nc) ;
	 for v = VAR
	    b_opt = zeros(1, X.nc) ;
	    Jv = strcmp(Y.vars, v) ; y = zscore(mean(Y.z(I,Jv), 2)) ;
	    [s_opt, b_opt(Jptr)] = crossvalidate(@lars, PAR.kfold, steps=100, zscore(X.x(I,Jptr)), y, PAR.reg) ;
	    Jeof = Jeof | b_opt != 0 ;
	 endfor

	 ## y = nansum([yP, yT, yP.*yT], 2) ;
	 ## [s_opt, b_opt(Jptr)] = crossvalidate(@lars, PAR.kfold, steps=100, zscore(X.x(I,Jptr)), y, PAR.reg) ;
	 ## Jeof = b_opt != 0 ;

      case "manual"
	 if !exist("NODISP", "var") || !ischar(NODISP), NODISP = "" ; endif 

	 NEOF = Y.nc:sum(Jptr) ;
	 if isempty(NEOF)
	    error("cal: too few predictors. %d < %d\n", sum(Jptr), Y.nc)
	 endif

	 [X.PJ C] = ptrSort(X.x(:,Jptr), X.J(Jptr), Y, PAR.ptrsel) ;
	 Xp = (1:X.nc)(Jptr)(X.PJ) ;
	 if PAR.dbg
	    Xj = Xjfun(X.J) ;
	    ptrlist = X.ptr(X.J(Xp)) ;
	    printf("\n   #:\tP ptr     \tno. \tcorr\t||  #:\tT ptr     \tno. \tcorr\n\n") ;
	    for j=1:length(ptrlist)/2
	       jP = 2*j-1 ; jT = 2*j ;
	       Pstr = sprintf("%-10s", ptrlist{jP}) ; Tstr = sprintf("%-10s", ptrlist{jT}) ;
	       XjP = Xj(Xp(jP)) ; XjT = Xj(Xp(jT)) ; Pc = abs(C(jP)) ; Tc = abs(C(jT)) ; 
	       printf("%4d:\t%s\t%3d\t%3.2f\t||%4d:\t%s\t%3d\t%3.2f\n", jP, Pstr, XjP, Pc, jT, Tstr, XjT, Tc) ;
	    endfor 
	 endif

	 JEOF = cell(1, length(NEOF)) ; JEOF = cellfun(@(j) Xp(1:j), num2cell(NEOF), "UniformOutput", false) ;
	 res = parcellfun_(PAR.ncpu, @(J) skl(J, Y, X), JEOF, PAR.par_opt{:}) ;
	 res = reshape(cell2mat(res), 2 + 2*PAR.kstest, [], length(NEOF)) ;

	 neof = neof_fig(NEOF, res, VAR) ;
	 Jeof = find2ind(Xp(1:neof), X.nc) ;

      case "stw"
	 if exist(file = [ldata "/ptr_stw.ob"])
	    load(file) ;
	 else
	    wX = X ; wX.x = X.x(:,Jptr) ;
	    stw = stw_fun(X, Y) ;
	    save(file, "stw") ;
	 endif
	 Jeof = stw{4} ;

      case "none"
	   Jeof = Jptr ;

      otherwise
	 warning("xds:xds", "cal: no ptr selection specified, using all predictors") ;
	 Jeof = true(1, X.nc) & Jptr ;

   endswitch

   if sum(Jeof) < Y.nc
      printf("%s:\t%d predictors selected, but at least %d (= number of predictands) are needed.\n", f(1).name, sum(Jeof), Y.nc) ;
      printf("%s:\tConsider using alternative PAR.optptr or PAR.trunc.\n", f(1).name) ;
      warning("xds:xds", "%s:\tContinuing with PAR.optptr = \"none\"", f(1).name) ;
      PAR.optptr = "none" ;
      Jeof = ptrsel (X, Y, Jptr) ;
   endif

   printf("\n%s: Selecting %d predictors\n", f(1).name, sum(Jeof)) ;
   printf("%s: --> %s\n", f(1).name, PAR.plist) ; 

endfunction


function [E, R] = do_cal (Y, X, Jeof, pcr_mode)

   ## usage: [E, R] = do_cal (Y, X, Jeof, pcr_mode)
   ## 
   ## calibrate models using Jeof predictors

   global PAR eof

   for is = 1:length(PAR.ssn)
      [Iy Ix] = selssn(PAR.ssn{is}, Y.id, X.id) ;
      wE = xrgr(Y.z(Iy,:), X.x(Ix,Jeof), "xds", Y.cov(:,:,is)) ;
      wR = xrgr(Y.z(Iy,:), X.x(Ix,Jeof), "rgr", Y.cov(:,:,is)) ;
      E(:,:,is) = R(:,:,is) = zeros(X.nc, Y.nc) ; E(Jeof,:,is) = wE ; R(Jeof,:,is) = wR ;
   endfor

   if nargin > 3 && pcr_mode
      E = mult(eof.E, E, 2, 1) ;
      R = mult(eof.E, R, 2, 1) ;
   endif

endfunction


function v = do_cov (u)

   ## usage: v = do_cov (u)
   ## 
   ## 

   global PAR

   v = u ;

   v.cal = "gregorian" ;

   for is = 1:length(PAR.ssn)

      I = sdate(v.id, PAR.cper(1,:), PAR.cper(2,:)) ;
      Iv = selssn(PAR.ssn{is}, v.id) ;
      c = cov(v.z(I&Iv,:)) ;

      if any(isnan(c(:)))
	 warning("xds:xds", "pdd: non-overlapping predictands. Guessing covariance.\n") ;
	 [i j] = find(isnan(c)) ;
	 n = length(i) ;
	 for k=1:n
	    I = 1:n != k ;
	    c(i(k),j(k)) = (mean(c(I,j(k)), 1) + mean(c(i(k),I), 2)) / 2 ;
	 endfor 
	 c = bend(c) ;
      endif 

      v.cov(:,:,is) = c ;

   endfor

endfunction
